# -*- coding: utf-8 -*-
"""
CIER 与 就业率（季度）的定性相关性分析（稳健读取版本）
输出：
  - cier_emp_quarterly.csv              # 季度对齐的 CIER、就业率及其 YoY
  - corr_cier_emp.csv                   # 水平 & YoY 的 Pearson/Spearman 相关与 p 值
  - overlay_levels.png / overlay_yoy.png
  - scatter_levels.png / scatter_yoy.png
  - rolling_corr_levels.png / rolling_corr_yoy.png   # 4 季滚动相关（定性）
  - ccf_levels.csv/png  ,  ccf_yoy.csv/png          # 跨相关（±6 季，k>0: CIER 领先）
  - seasonal_quarter_stats.csv                      # CIER 与就业率的季度季节性
  - qual_summary.csv                                 # 符号一致率、最长同向/反向共振区间
"""

import os
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import pearsonr, spearmanr

# ========= 改这里 =========
PATH_CIER = "/Users/linshangjin/25CCM/NKU-C/dataset_csv/2016-2025 CIER.csv"       # 你的 CIER 季度数据（CSV）
PATH_EMP  = "/Users/linshangjin/25CCM/NKU-C/dataset_csv/data_processed.xlsx"      # 就业率日/半月数据（Excel，多 sheet）
OUT_DIR   = "/Users/linshangjin/25CCM/NKU-C/t2/CIER"       # 输出目录
# =========================

os.makedirs(OUT_DIR, exist_ok=True)

# ---------- 读取工具（稳健） ----------
def read_csv_flexible(path):
    # 优先 utf-8，失败再 gbk
    try:
        return pd.read_csv(path, encoding="utf-8")
    except UnicodeDecodeError:
        return pd.read_csv(path, encoding="gbk")

def detect_date_column(cols):
    cand = ["date", "时间", "period", "季度", "Quarter", "quarter", "Unnamed: 0"]
    for c in cand:
        if c in cols:
            return c
    # 兜底：找包含 'date' 或 'time' 的
    for c in cols:
        lc = str(c).lower()
        if "date" in lc or "time" in lc or "quarter" in lc:
            return c
    return None

def detect_value_column(df, preferred=("CIER","CIER指数","value","指数","val")):
    for c in preferred:
        if c in df.columns:
            return c
    # 找到首个“数值型且非日期”的列
    for c in df.columns:
        if pd.api.types.is_numeric_dtype(df[c]):
            return c
    # 尝试把非数值列转一下
    for c in df.columns:
        if c != detect_date_column(df.columns):
            try:
                df[c] = pd.to_numeric(df[c], errors="coerce")
                if df[c].notna().sum() > 0:
                    return c
            except Exception:
                continue
    return None

def read_cier_quarterly(path):
    df = read_csv_flexible(path)
    dcol = detect_date_column(df.columns)
    if dcol is None:
        raise ValueError("CIER 文件未找到日期列（如 date/时间/季度/period/Unnamed: 0）")
    df[dcol] = pd.to_datetime(df[dcol])
    vcol = detect_value_column(df)
    if vcol is None:
        raise ValueError("CIER 文件未找到数值列（如 CIER/CIER指数/value）")
    s = df[[dcol, vcol]].dropna().sort_values(dcol).set_index(dcol)[vcol].astype(float)
    # 规范为季度：取每季最后一个观察
    s_q = s.resample("Q").last().rename("CIER")
    return s_q

def read_employment_quarterly(path):
    xl = pd.ExcelFile(path)
    cand_val = ["employment_rate","就业率","y_pred","y","employment","rate"]
    parts = []
    found = False
    for sh in xl.sheet_names:
        t = xl.parse(sh)
        # 找日期列
        dcol = None
        for c in ["date","日期","time","Time","时间"]:
            if c in t.columns:
                dcol = c; break
        if dcol is None: 
            continue
        # 找就业率列
        vcol = None
        for c in cand_val:
            if c in t.columns:
                vcol = c; break
        if vcol is None:
            # 兜底：找非日期的第一个数值列
            for c in t.columns:
                if c != dcol and pd.api.types.is_numeric_dtype(t[c]):
                    vcol = c; break
        if vcol is None: 
            continue
        tmp = t[[dcol, vcol]].copy()
        tmp[dcol] = pd.to_datetime(tmp[dcol])
        tmp[vcol] = pd.to_numeric(tmp[vcol], errors="coerce")
        parts.append(tmp.rename(columns={dcol:"date", vcol:"employment_rate"}))
        found = True
    if not found:
        raise ValueError("未在 data_processed.xlsx 中找到包含日期和就业率的工作表")
    emp = pd.concat(parts, ignore_index=True).dropna(subset=["date"]).sort_values("date")
    emp = emp.set_index("date")
    # 聚合到季度：均值
    y_q = emp["employment_rate"].resample("Q").mean().rename("EMP")
    return y_q

def yoy_quarter(s, k=4, log=False):
    if log:
        s = s.copy()
        if (s <= 0).any():
            return s.diff(k)  # 非正值不适合对数
        return np.log(s).diff(k)
    return s.diff(k)

def ccf_xy(x, y, maxlag=6):
    """跨相关：k>0 表示 x 领先 y"""
    x = (x - x.mean())/x.std()
    y = (y - y.mean())/y.std()
    out = {}
    for k in range(-maxlag, maxlag+1):
        if k >= 0:
            xs, ys = x.iloc[:-k or None], y.iloc[k:] if k else y
        else:
            xs, ys = x.iloc[-k:], y.iloc[:k or None]
        if len(xs) > 3:
            out[k] = float(np.corrcoef(xs, ys)[0,1])
    return out

# ---------- 读取与对齐 ----------
cier = read_cier_quarterly(PATH_CIER)
emp_q = read_employment_quarterly(PATH_EMP)
df = pd.concat([cier, emp_q], axis=1).dropna().rename(columns={"CIER":"CIER","EMP":"EMP"})
# YoY（季度同比差分 4 期；指数/比率都用差分以避免伪高相关）
df["CIER_YoY"] = yoy_quarter(df["CIER"], k=4, log=False)
df["EMP_YoY"]  = yoy_quarter(df["EMP"],  k=4, log=False)
df.to_csv(os.path.join(OUT_DIR, "cier_emp_quarterly.csv"))

# ---------- 相关（水平 & YoY） ----------
def corr_pair(a, b):
    idx = a.dropna().index.intersection(b.dropna().index)
    if len(idx) < 6:
        return np.nan, np.nan, np.nan, np.nan, int(len(idx))
    rP, pP = pearsonr(a.loc[idx], b.loc[idx])
    rS, pS = spearmanr(a.loc[idx], b.loc[idx])
    return float(rP), float(pP), float(rS), float(pS), int(len(idx))

rows = []
rP,pP,rS,pS,n = corr_pair(df["CIER"], df["EMP"])
rows.append({"type":"levels","pearson":rP,"pearson_p":pP,"spearman":rS,"spearman_p":pS,"n_obs":n})
rP,pP,rS,pS,n = corr_pair(df["CIER_YoY"], df["EMP_YoY"])
rows.append({"type":"yoy(Δ4)","pearson":rP,"pearson_p":pP,"spearman":rS,"spearman_p":pS,"n_obs":n})
pd.DataFrame(rows).to_csv(os.path.join(OUT_DIR, "corr_cier_emp.csv"), index=False)

# ---------- 叠线图（水平 & YoY） ----------
fig = plt.figure()
df[["CIER","EMP"]].plot(ax=plt.gca())
plt.title("CIER vs Employment (Quarterly Levels)")
plt.xlabel("Quarter"); plt.ylabel("Level / Index")
plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, "overlay_levels.png"), dpi=150); plt.close(fig)

fig = plt.figure()
df[["CIER_YoY","EMP_YoY"]].dropna().plot(ax=plt.gca())
plt.axhline(0, linewidth=1)
plt.title("CIER vs Employment (Quarterly YoY, Δ4)")
plt.xlabel("Quarter"); plt.ylabel("YoY (Δ4)")
plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, "overlay_yoy.png"), dpi=150); plt.close(fig)

# ---------- 散点图（水平 & YoY） ----------
def scatter_with_fit(x, y, xlabel, ylabel, path):
    fig = plt.figure()
    xx, yy = x.values, y.values
    plt.scatter(xx, yy, s=30)
    if len(xx) >= 3:
        a, b = np.polyfit(xx, yy, 1)
        grid = np.linspace(np.nanmin(xx), np.nanmax(xx), 100)
        plt.plot(grid, a*grid + b)
    plt.xlabel(xlabel); plt.ylabel(ylabel)
    plt.tight_layout(); plt.savefig(path, dpi=150); plt.close(fig)

pairL = df[["CIER","EMP"]].dropna()
if len(pairL) >= 3:
    scatter_with_fit(pairL["CIER"], pairL["EMP"], "CIER (level)", "EMP (level)",
                     os.path.join(OUT_DIR, "scatter_levels.png"))

pairY = df[["CIER_YoY","EMP_YoY"]].dropna()
if len(pairY) >= 3:
    scatter_with_fit(pairY["CIER_YoY"], pairY["EMP_YoY"], "CIER YoY (Δ4)", "EMP YoY (Δ4)",
                     os.path.join(OUT_DIR, "scatter_yoy.png"))

# ---------- 4 季滚动相关（定性） ----------
def rolling_corr(a, b, w=4):
    return a.rolling(w).corr(b)

rollL = rolling_corr(df["CIER"], df["EMP"], w=4)
rollY = rolling_corr(df["CIER_YoY"], df["EMP_YoY"], w=4)

fig = plt.figure()
rollL.plot(ax=plt.gca())
plt.title("Rolling Corr (w=4): Levels")
plt.xlabel("Quarter"); plt.ylabel("Corr")
plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, "rolling_corr_levels.png"), dpi=150); plt.close(fig)

fig = plt.figure()
rollY.plot(ax=plt.gca())
plt.title("Rolling Corr (w=4): YoY (Δ4)")
plt.xlabel("Quarter"); plt.ylabel("Corr")
plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, "rolling_corr_yoy.png"), dpi=150); plt.close(fig)

# ---------- CCF（±6 季；k>0 表示 CIER 领先） ----------
def save_ccf(x, y, fname_csv, fname_png, title):
    c = ccf_xy(x.dropna(), y.dropna(), maxlag=6)
    cdf = pd.DataFrame({"lag_q": list(c.keys()), "corr": list(c.values())}).sort_values("lag_q")
    cdf.to_csv(os.path.join(OUT_DIR, fname_csv), index=False)
    fig = plt.figure()
    plt.bar(cdf["lag_q"], cdf["corr"])
    plt.title(title)
    plt.xlabel("Lag (quarters)   k>0: CIER leads"); plt.ylabel("Correlation")
    plt.tight_layout(); plt.savefig(os.path.join(OUT_DIR, fname_png), dpi=150); plt.close(fig)

if len(pairL) >= 8:
    save_ccf(df["CIER"], df["EMP"], "ccf_levels.csv", "ccf_levels.png",
             "CCF (Levels): CIER vs EMP")

if len(pairY) >= 8:
    save_ccf(df["CIER_YoY"], df["EMP_YoY"], "ccf_yoy.csv", "ccf_yoy.png",
             "CCF (YoY Δ4): CIER vs EMP")

# ---------- 季节性统计（Q1~Q4） ----------
season = pd.DataFrame({
    "CIER": df["CIER"],
    "EMP": df["EMP"],
    "Q": df.index.quarter
}).dropna()

def q_stats(s):
    return s.groupby(season["Q"]).agg(mean="mean", median="median",
                                      p25=lambda x: x.quantile(0.25),
                                      p75=lambda x: x.quantile(0.75),
                                      count="count")

stats = pd.concat({
    "CIER": q_stats(season["CIER"]),
    "EMP":  q_stats(season["EMP"])
}, axis=1)
stats.to_csv(os.path.join(OUT_DIR, "seasonal_quarter_stats.csv"))

# ---------- 定性摘要：符号一致率 + 最长共振区间 ----------
qual_rows = []
def sign_co_move(x, y, name):
    idx = x.dropna().index.intersection(y.dropna().index)
    x1, y1 = x.loc[idx], y.loc[idx]
    if len(idx) == 0:
        return {"variable": name, "sign_agree": np.nan}
    agree = float((np.sign(x1) == np.sign(y1)).mean())
    # 找最长连续同向（>0）或同向（<0）的片段
    state = np.where((x1>0)&(y1>0), 1, np.where((x1<0)&(y1<0), -1, 0))
    segs = []; start=None; cur=0
    for i, s in enumerate(state):
        if s != 0 and start is None:
            start = i; cur=s
        elif s != 0 and s != cur:
            segs.append((start, i-1, cur)); start=i; cur=s
        elif s == 0 and start is not None:
            segs.append((start, i-1, cur)); start=None
    if start is not None:
        segs.append((start, len(state)-1, cur))
    segs = sorted(segs, key=lambda t: (t[1]-t[0]+1), reverse=True)
    rec = {"variable": name, "sign_agree": agree}
    if segs:
        a,b,sgn = segs[0]
        idxv = list(idx)
        rec.update({
            "longest_seg_start": str(idxv[a].date()),
            "longest_seg_end":   str(idxv[b].date()),
            "longest_seg_len_q": int(b-a+1),
            "longest_seg_type":  "both_up" if sgn==1 else "both_down"
        })
    return rec

qual_rows.append(sign_co_move(df["CIER"], df["EMP"], "levels"))
qual_rows.append(sign_co_move(df["CIER_YoY"], df["EMP_YoY"], "yoy(Δ4)"))
pd.DataFrame(qual_rows).to_csv(os.path.join(OUT_DIR, "qual_summary.csv"), index=False)

print("✅ Done. 输出目录：", os.path.abspath(OUT_DIR))
